package com.example.testapp;

import android.util.Log;

public class Solver {
	private int N;
	private int K;
	private Functions data;
	private double h;
	
	static final int exp_1 = 0;
	static final int exp_2_1 = 1;
	static final int exp_3_2 = 2;
	static final int exp_2_2 = 3;
	
	static final int imp_1 = 4;
	static final int imp_2_1 = 5;
	static final int imp_2_2 = 6;
	static final int imp_3_2 = 7;
	static final int real = 8;
	
	public Solver(int _N, int _K)
	{
		N = _N; K = _K;
	}
	
	public double getH()
	{
		return h;
	}
	
	public void setN(int _N)
	{
		N = _N;
	}
	
	public void setK(int _K)
	{
		K = _K;
	}
	
	public int getN()
	{
		return N;
	}
	
	public int getK()
	{
		return K;
	}
	
	public double[][] solve(int type, Functions _data, double T)
	{
		data = _data;
		switch(type)
		{
		case exp_1:
			return explicit_1(T);
		case exp_2_1:
			return explicit_2_1(T);
		case exp_3_2:
			return explicit_3_2(T);
		case exp_2_2:
			return explicit_2_2(T);
		case imp_1:
			return implicit_1(T);
		case imp_2_1:
			return implicit_2_1(T);
		case imp_2_2:
			return implicit_2_2(T);
		case imp_3_2:
			return implicit_3_2(T);
		}
		
		return new double[1][1];
	}
	
	private double [] sweep(int n, double a[], double c[], double b[], double f[])
	{
		double x[] = new double[n];
		double m;
		for(int i = 1; i < n; ++i)
		{
			m = a[i]/c[i-1];
			c[i] = c[i] - m*b[i-1];
			f[i] = f[i] - m*f[i-1];
		}
		x[n-1] = f[n-1]/c[n-1];
		for(int i = n-2; i >= 0; --i)
		{
			x[i] = (f[i]-b[i]*x[i+1])/c[i];
		}
		return x;
	}
	
	public double[][] explicit_1(double T)
	{
		double[][] U = new double[N+1][K+1];
		for(int i = 0; i <= N; ++i)
			for(int k = 0; k <= K; ++k)
				U[i][k] = 0;
		double tau = T/K;
		h = data.L/N;
		
		double sigma = (data.A*data.A*tau*tau)/(h*h);
		
        Log.d("exp", "sigma" + Double.toString(sigma) +
        		" tau " + Double.toString(tau) + " K " +
        		Integer.toString(K));

//	    while(sigma > 0.5)
//	    {
//	        tau -= 0.000001;
//	        K=(int)(T/tau);
//	        sigma = (data.A*data.A*tau*tau)/(h*h);
//	        Log.d("exp", "sigma" + Double.toString(sigma) +
//	        		" tau " + Double.toString(tau) + " K " +
//	        		Integer.toString(K));
//	    }
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2.;
		}
		for(int k = 0; k <=K ; ++k)
		{
			U[0][k] = data.phi_0(data.t0 + tau * k);
			U[N][k] = data.phi_l(data.t0 + tau * k);
		}
		for(int k = 1; k < K; ++k)
			for(int i = 1; i < N; ++i)
			{
				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2.*U[i][k] + U[i-1][k])/(h * h) +
						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2. * h) +
						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i) + 
						2. * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
						data.D * U[i][k]/tau)/(1./(tau * tau) + data.D/tau);
//				U[i][k+1] = ((data.A * data.A * tau * tau)*(U[i+1][k] - 2.*U[i][k] + U[i-1][k])/(h * h) +
//						(data.B * data.B * tau * tau)*(U[i+1][k] - U[i-1][k])/(2. * h) +
//						data.C * tau * tau * U[i][k] + tau * tau * data.f(data.x0 + h * i, data.t0 + tau * i)) /
//						(1. + tau) + 
//						2. * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
//						data.D * U[i][k]/tau;
//				U[i][k+1] = U[i+1][k]*(data.A*data.A*tau*tau/(h*h) + data.B*tau*tau/(2*h)) + U[i][k]*(-2*data.A*data.A*tau*tau/(h*h) + 2 +
//						data.C*tau*tau) + U[i-1][k]*(data.A*data.A*tau*tau/(h*h) - data.B*tau*tau/(2*h)) - U[i][k-1] + tau*tau*data.f(data.x0 + i*h, data.t0 + k*tau);
				//if(Double.isNaN(U[i][k+1])) break;
				//Log.d("explcit_1", Double.toString(U[i][k+1]));
				//Log.d("numbers", Integer.toString(i) + " " + Integer.toString(k+1));
			}
		return U;
	}
	
	public void setData(Functions _data)
	{
		data = _data;
	}
	
	public double[][] explicit_2_1(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}

		for(int k = 1; k < K; ++k)
		{
			for(int i = 1; i < N; ++i)
			{
//				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2*U[i][k] + U[i-1][k])/(h * h) +
//						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2 * h) +
//						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i)) /
//						(1/(tau * tau) + data.D/tau) + 
//						2 * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
//						data.D * U[i][k]/tau;
//				U[i][k+1] = U[i+1][k]*(data.A*data.A*tau*tau/(h*h) + data.B*tau*tau/(2*h)) + U[i][k]*(-2*data.A*data.A*tau*tau/(h*h) + 2 +
//						data.C*tau*tau) + U[i-1][k]*(data.A*data.A*tau*tau/(h*h) - data.B*tau*tau/(2*h)) - U[i][k-1] + tau*tau*data.f(data.x0 + i*h, data.t0 + k*tau);
				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2.*U[i][k] + U[i-1][k])/(h * h) +
						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2. * h) +
						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i) + 
						2. * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
						data.D * U[i][k]/tau)/(1./(tau * tau) + data.D/tau);
			}
			U[0][k+1] = -(data.Alpha/h) * U[1][k+1] / (data.Beta - data.Alpha/h) + 
					data.phi_0(data.t0 + tau * (k+1)) / (data.Beta - data.Alpha/h);
			U[N][k+1] = (data.Gamma/h)*U[N-1][k+1] / (data.Delta + data.Gamma/h) +
					data.phi_l(data.t0 + tau * (k+1))/(data.Delta + data.Gamma/h);
		}
		
		return U;
	}
	
	public double[][] explicit_3_2(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}

		for(int k = 1; k < K; ++k)
		{
			for(int i = 1; i < N; ++i)
			{
//				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2*U[i][k] + U[i-1][k])/(h * h) +
//						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2 * h) +
//						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i)) /
//						(1/(tau * tau) + data.D/tau) + 
//						2 * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
//						data.D * U[i][k]/tau;
//				U[i][k+1] = U[i+1][k]*(data.A*data.A*tau*tau/(h*h) + data.B*tau*tau/(2*h)) + U[i][k]*(-2*data.A*data.A*tau*tau/(h*h) + 2 +
//						data.C*tau*tau) + U[i-1][k]*(data.A*data.A*tau*tau/(h*h) - data.B*tau*tau/(2*h)) - U[i][k-1] + tau*tau*data.f(data.x0 + i*h, data.t0 + k*tau);
				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2.*U[i][k] + U[i-1][k])/(h * h) +
						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2. * h) +
						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i) + 
						2. * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
						data.D * U[i][k]/tau)/(1./(tau * tau) + data.D/tau);
			}
		U[0][k+1] = (data.phi_0(data.t0 + tau * (k+1)) - data.Alpha*(4*U[1][k+1] - U[2][k+1])/(2*h)) /
				(-3*data.Alpha/(2*h) + data.Beta);
		U[N][k+1] = ( data.phi_l(data.t0 + tau * (k+1)) - data.Gamma*(U[N-2][k+1] - 4*U[N-1][k+1])/(2*h) ) / 
				(3*data.Gamma/(2*h) + data.Delta);
		}
		
		return U;
	}
	
	public double[][] explicit_2_2(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}

		for(int k = 1; k < K; ++k)
		{
			for(int i = 1; i < N; ++i)
			{
//				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2*U[i][k] + U[i-1][k])/(h * h) +
//						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2 * h) +
//						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i)) /
//						(1/(tau * tau) + data.D/tau) + 
//						2 * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
//						data.D * U[i][k]/tau;
//				U[i][k+1] = U[i+1][k]*(data.A*data.A*tau*tau/(h*h) + data.B*tau*tau/(2*h)) + U[i][k]*(-2*data.A*data.A*tau*tau/(h*h) + 2 +
//						data.C*tau*tau) + U[i-1][k]*(data.A*data.A*tau*tau/(h*h) - data.B*tau*tau/(2*h)) - U[i][k-1] + tau*tau*data.f(data.x0 + i*h, data.t0 + k*tau);
				U[i][k+1] = ((data.A * data.A)*(U[i+1][k] - 2.*U[i][k] + U[i-1][k])/(h * h) +
						(data.B * data.B)*(U[i+1][k] - U[i-1][k])/(2. * h) +
						data.C * U[i][k] + data.f(data.x0 + h * i, data.t0 + tau * i) + 
						2. * U[i][k]/(tau * tau) - U[i][k-1]/(tau * tau) + 
						data.D * U[i][k]/tau)/(1./(tau * tau) + data.D/tau);
			}
		
			U[0][k+1] = ( data.phi_0(data.t0 + tau * (k+1)) - data.Alpha/(2*data.A*data.A-data.B*h) * 
					(2*data.A*data.A/h * U[1][k+1] + 2*h/(tau*tau)*U[0][k] - 
							h/(tau*tau)*U[0][k-1] + data.D*h/tau * U[0][k] +
							h*data.f(data.x0, data.t0 + tau * (k+1)) ) ) / 
							( data.Alpha/(2*data.A*data.A - data.B*h) * 
							( -2*data.A*data.A/h - h/(tau*tau) - data.D*h/tau + data.C*h ) + data.Beta );
			U[N][k+1] = ( data.phi_l(data.t0 + tau * (k+1)) - data.Gamma/(2*data.A*data.A+data.B*h) * 
					(-2*data.A*data.A/h * U[N-1][k+1] + h/(tau*tau)*(U[N][k-1] - 
							2*U[N][k]) - data.D*h/tau * U[N][k] -
							h*data.f(data.L, data.t0 + tau * (k+1)) ) ) / 
							( data.Gamma/(2*data.A*data.A + data.B*h) * 
							( 2*data.A*data.A/h + h/(tau*tau) + data.D*h/tau - data.C*h ) + data.Delta );
		}
		
		return U;
	}
	
	public double[][] implicit_1(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}
		
		for(int k = 2; k <=K ; ++k)
		{
			U[0][k] = data.phi_0(data.t0 + tau * k);
			U[N][k] = data.phi_l(data.t0 + tau * k);
		}
		
		for(int k = 1; k < K; ++k)
		{
			double A[] = new double[N-1];
			double B[] = new double[N-1];
			double C[] = new double[N-1];
			double F[] = new double[N-1];
			
			A[0] = 0;
			C[0] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C);
			B[0] = data.A*data.A/(h*h) + data.B/(2*h);
			F[0] = 1/(tau*tau)*U[1][k-1] - (2/(tau*tau) + data.D/tau)*U[1][k] - data.f(data.x0 + h, data.t0 + tau) - 
					(data.A*data.A/(h*h) - data.B/(2*h)) * data.phi_0(data.t0 + (k+1) * tau);
			for(int j = 2; j < N-1; ++j)
			{
				A[j-1] = data.A*data.A/(h*h) - data.B/(2*h);
				C[j-1] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C);
				B[j-1] = data.A*data.A/(h*h) + data.B/(2*h);
				F[j-1] = 1/(tau*tau)*U[j][k-1] - (2/(tau*tau) + data.D/tau)*U[j][k] - data.f(data.x0 + j*h, data.t0 + (k+1)*tau);
			}
			
			A[N-2] = data.A*data.A/(h*h) - data.B/(2*h);
			C[N-2] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C);
			B[N-2] = 0;
			F[N-2] = 1/(tau*tau)*U[N-1][k-1] - (2/(tau*tau) + data.D/tau)*U[N-1][k] - data.f(data.x0 + (N-1)*h, data.t0 + (k+1)*tau) -
					(data.A * data.A/(h*h) + data.B/(2*h))*data.phi_l(data.t0 + (k+1) * tau);
			
			double X[] = sweep(N-1, A, C, B, F);
			for(int i= 0; i < N-1; ++i)
			{
				U[i+1][k+1] = X[i];
				//Log.d("implicit_1", Double.toString(U[i+1][k+1]));
				//Log.d("numbers", Integer.toString(i) + " " + Integer.toString(k+1));
			}
		}
		
		return U;
	}
	
	public double[][] implicit_2_1(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}
		
		for(int k = 1; k < K; ++k)
		{
			double A[] = new double[N-1];
			double B[] = new double[N-1];
			double C[] = new double[N-1];
			double F[] = new double[N-1];
			
			A[0] = 0;
			C[0] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C - 
					(data.A*data.A/(h*h) - data.B/(2*h))*((data.Alpha/h)/(data.Beta - 
							data.Alpha/h) ) );
			B[0] = data.A*data.A/(h*h) + data.B/(2*h);
			F[0] = 1/(tau*tau)*U[1][k-1] - (2/(tau*tau) + data.D/tau)*U[1][k] - data.f(data.x0 + h, data.t0 + tau) - 
					(data.A*data.A/(h*h) - data.B/(2*h)) * (
					///////////////////////////////////// 
					data.phi_0(data.t0 + tau * (k+1)) / (data.Beta - data.Alpha/h) 
					);
			////////////////////////////////////////////////////
					
			for(int j = 2; j < N-1; ++j)
			{
				A[j-1] = data.A*data.A/(h*h) - data.B/(2*h);
				C[j-1] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C);
				B[j-1] = data.A*data.A/(h*h) + data.B/(2*h);
				F[j-1] = 1/(tau*tau)*U[j][k-1] - (2/(tau*tau) + data.D/tau)*U[j][k] - data.f(data.x0 + j*h, data.t0 + (k+1)*tau);
			}
			
			A[N-2] = data.A*data.A/(h*h) - data.B/(2*h) ;
			C[N-2] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C - 
					(data.A*data.A/(h*h) + data.B/(2*h) ) *
					( (data.Gamma/h)/(data.Delta + data.Gamma/h) ));
			B[N-2] = 0;
			F[N-2] = 1/(tau*tau)*U[N-1][k-1] - (2/(tau*tau) + data.D/tau)*U[N-1][k] - data.f(data.x0 + (N-1)*h, data.t0 + (k+1)*tau) -
					(data.A * data.A/(h*h) + data.B/(2*h))* (
			///////////////////////////////////////////////////////
					data.phi_l(data.t0 + tau * (k+1))/(data.Delta + data.Gamma/h)
			//////////////////////////////////////////////////////////////////
					);
			
			double X[] = sweep(N-1, A, C, B, F);
			for(int i= 0; i < N-1; ++i)
			{
				U[i+1][k+1] = X[i];
			}
			
			U[0][k+1] = -(data.Alpha/h) * U[1][k+1] / (data.Beta - data.Alpha/h) + 
					data.phi_0(data.t0 + tau * (k+1)) / (data.Beta - data.Alpha/h);
			U[N][k+1] = (data.Gamma/h)*U[N-1][k+1] / (data.Delta + data.Gamma/h) +
					data.phi_l(data.t0 + tau * (k+1))/(data.Delta + data.Gamma/h);
		}
		
		return U;
	}
	
	public double[][] implicit_2_2(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}
		
		for(int k = 1; k < K; ++k)
		{
			double A[] = new double[N-1];
			double B[] = new double[N-1];
			double C[] = new double[N-1];
			double F[] = new double[N-1];
			
			A[0] = 0;
			C[0] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C 
					- ( (data.A*data.A/(h*h) - data.B/(2*h)) * 
					( data.Alpha*2*data.A*data.A/(h*(2*data.A*data.A - data.B*h)) ) /
					( data.Alpha/(2*data.A*data.A - data.B*h) * 
							( -2*data.A*data.A/h - h/(tau*tau) - data.D*h/tau + data.C*h ) + data.Beta )
							));
			B[0] = data.A*data.A/(h*h) + data.B/(2*h);
			F[0] = 1/(tau*tau)*U[1][k-1] - (2/(tau*tau) + data.D/tau)*U[1][k] - data.f(data.x0 + h, data.t0 + tau) - 
						(data.A*data.A/(h*h) - data.B/(2*h)) * (
						//////////////////////////////////////////////////////////
						( data.phi_0(data.t0 + tau * (k+1)) - data.Alpha/(2*data.A*data.A-data.B*h) * 
						( 2*h/(tau*tau)*U[0][k] - 
						h/(tau*tau)*U[0][k-1] + data.D*h/tau * U[0][k] +
						h*data.f(data.x0, data.t0 + tau * (k+1)) ) ) / 
						( data.Alpha/(2*data.A*data.A - data.B*h) * 
						( -2*data.A*data.A/h - h/(tau*tau) - data.D*h/tau + data.C*h ) + data.Beta )
							////////////////////////////////////////////////////////
							);
			for(int j = 2; j < N-1; ++j)
			{
				A[j-1] = data.A*data.A/(h*h) - data.B/(2*h);
				C[j-1] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C);
				B[j-1] = data.A*data.A/(h*h) + data.B/(2*h);
				F[j-1] = 1/(tau*tau)*U[j][k-1] - (2/(tau*tau) + data.D/tau)*U[j][k] - data.f(data.x0 + j*h, data.t0 + (k+1)*tau);
			}
			
			A[N-2] = data.A*data.A/(h*h) - data.B/(2*h);
			C[N-2] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C +
					(data.A * data.A/(h*h) + data.B/(2*h))*( data.Gamma/(2*data.A*data.A+data.B*h)*
					2*data.A*data.A/h) / ( data.Gamma/(2*data.A*data.A + data.B*h) * 
							( 2*data.A*data.A/h + h/(tau*tau) + data.D*h/tau - data.C*h ) + data.Delta ));
			B[N-2] = 0;
			F[N-2] = 1/(tau*tau)*U[N-1][k-1] - (2/(tau*tau) + data.D/tau)*U[N-1][k] - data.f(data.x0 + (N-1)*h, data.t0 + (k+1)*tau) -
					(data.A * data.A/(h*h) + data.B/(2*h))* (
							//////////////////////////////////////////////////////////////
							( data.phi_l(data.t0 + tau * (k+1)) - data.Gamma/(2*data.A*data.A+data.B*h) * 
							( h/(tau*tau)*(U[N][k-1] - 
							2*U[N][k]) - data.D*h/tau * U[N][k] -
							h*data.f(data.L, data.t0 + tau * (k+1)) ) ) / 
							( data.Gamma/(2*data.A*data.A + data.B*h) * 
							( 2*data.A*data.A/h + h/(tau*tau) + data.D*h/tau - data.C*h ) + data.Delta )
							//////////////////////////////////////////////////////////////////
							);
			
			double X[] = sweep(N-1, A, C, B, F);
			for(int i= 0; i < N-1; ++i)
			{
				U[i+1][k+1] = X[i];
			}
			
			U[0][k+1] = ( data.phi_0(data.t0 + tau * (k+1)) - data.Alpha/(2*data.A*data.A-data.B*h) * 
					(2*data.A*data.A/h * U[1][k+1] + 2*h/(tau*tau)*U[0][k] - 
							h/(tau*tau)*U[0][k-1] + data.D*h/tau * U[0][k] +
							h*data.f(data.x0, data.t0 + tau * (k+1)) ) ) / 
							( data.Alpha/(2*data.A*data.A - data.B*h) * 
							( -2*data.A*data.A/h - h/(tau*tau) - data.D*h/tau + data.C*h ) + data.Beta );
			U[N][k+1] = ( data.phi_l(data.t0 + tau * (k+1)) - data.Gamma/(2*data.A*data.A+data.B*h) * 
					(-2*data.A*data.A/h * U[N-1][k+1] + h/(tau*tau)*(U[N][k-1] - 
							2*U[N][k]) - data.D*h/tau * U[N][k] -
							h*data.f(data.L, data.t0 + tau * (k+1)) ) ) / 
							( data.Gamma/(2*data.A*data.A + data.B*h) * 
							( 2*data.A*data.A/h + h/(tau*tau) + data.D*h/tau - data.C*h ) + data.Delta );
		}
		
		return U;
	}
	
	public double[][] implicit_3_2(double T)
	{
		double[][] U = new double[N+1][K+1];
		double tau = T/K;
		h = data.L/N;
		
		for(int i = 0; i <= N; ++i)
		{
			U[i][0] = data.psy_1(data.x0 + h * i);
			U[i][1] = data.psy_1(data.x0 + h * i) + data.psy_2(data.x0 + h * i) * tau +
					(data.A * data.A * data.diff2_psy_1(data.x0 + h * i) +
							data.B * data.diff1_psy_1(data.x0 + h * i) + 
							data.C * data.psy_1(data.x0 + h * i) + 
							data.f(data.x0 + h * i, data.t0 + tau) - 
							data.D * data.diff1_psy_2(data.x0 + h * i)) * 
							(tau * tau)/2;
		}
		
		for(int k = 1; k < K; ++k)
		{
			double A[] = new double[N-1];
			double B[] = new double[N-1];
			double C[] = new double[N-1];
			double F[] = new double[N-1];
			
			A[0] = 0;
			C[0] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C 
					- ( (data.A*data.A/(h*h) - data.B/(2*h)) * 
					( data.Alpha*2*data.A*data.A/(h*(2*data.A*data.A - data.B*h)) ) /
					( data.Alpha/(2*data.A*data.A - data.B*h) * 
							( -2*data.A*data.A/h - h/(tau*tau) - data.D*h/tau + data.C*h ) + data.Beta )
							));
			B[0] = data.A*data.A/(h*h) + data.B/(2*h);
			F[0] = 1/(tau*tau)*U[1][k-1] - (2/(tau*tau) + data.D/tau)*U[1][k] - data.f(data.x0 + h, data.t0 + tau) - 
						(data.A*data.A/(h*h) - data.B/(2*h)) * (
						//////////////////////////////////////////////////////////
						( data.phi_0(data.t0 + tau * (k+1)) - data.Alpha/(2*data.A*data.A-data.B*h) * 
						( 2*h/(tau*tau)*U[0][k] - 
						h/(tau*tau)*U[0][k-1] + data.D*h/tau * U[0][k] +
						h*data.f(data.x0, data.t0 + tau * (k+1)) ) ) / 
						( data.Alpha/(2*data.A*data.A - data.B*h) * 
						( -2*data.A*data.A/h - h/(tau*tau) - data.D*h/tau + data.C*h ) + data.Beta )
							////////////////////////////////////////////////////////
							);
			for(int j = 2; j < N-1; ++j)
			{
				A[j-1] = data.A*data.A/(h*h) - data.B/(2*h);
				C[j-1] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C);
				B[j-1] = data.A*data.A/(h*h) + data.B/(2*h);
				F[j-1] = 1/(tau*tau)*U[j][k-1] - (2/(tau*tau) + data.D/tau)*U[j][k] - data.f(data.x0 + j*h, data.t0 + (k+1)*tau);
			}
			
			A[N-2] = data.A*data.A/(h*h) - data.B/(2*h);
			C[N-2] = -(1/(tau*tau) + data.D/tau + 2*data.A*data.A/(h*h) - data.C +
					(data.A * data.A/(h*h) + data.B/(2*h))*( data.Gamma/(2*data.A*data.A+data.B*h)*
					2*data.A*data.A/h) / ( data.Gamma/(2*data.A*data.A + data.B*h) * 
							( 2*data.A*data.A/h + h/(tau*tau) + data.D*h/tau - data.C*h ) + data.Delta ));
			B[N-2] = 0;
			F[N-2] = 1/(tau*tau)*U[N-1][k-1] - (2/(tau*tau) + data.D/tau)*U[N-1][k] - data.f(data.x0 + (N-1)*h, data.t0 + (k+1)*tau) -
					(data.A * data.A/(h*h) + data.B/(2*h))* (
							//////////////////////////////////////////////////////////////
							( data.phi_l(data.t0 + tau * (k+1)) - data.Gamma/(2*data.A*data.A+data.B*h) * 
							( h/(tau*tau)*(U[N][k-1] - 
							2*U[N][k]) - data.D*h/tau * U[N][k] -
							h*data.f(data.L, data.t0 + tau * (k+1)) ) ) / 
							( data.Gamma/(2*data.A*data.A + data.B*h) * 
							( 2*data.A*data.A/h + h/(tau*tau) + data.D*h/tau - data.C*h ) + data.Delta )
							//////////////////////////////////////////////////////////////////
							);
			
			double X[] = sweep(N-1, A, C, B, F);
			for(int i= 0; i < N-1; ++i)
			{
				U[i+1][k+1] = X[i];
			}
			
			U[0][k+1] = (data.phi_0(data.t0 + tau * (k+1)) - data.Alpha*(4*U[1][k+1] - U[2][k+1])/(2*h)) /
					(-3*data.Alpha/(2*h) + data.Beta);
			U[N][k+1] = ( data.phi_l(data.t0 + tau * (k+1)) - data.Gamma*(U[N-2][k+1] - 4*U[N-1][k+1])/(2*h) ) / 
					(3*data.Gamma/(2*h) + data.Delta);
		}
		
		return U;
	}
	
	public double[] uReal(double T)
	{
		h = data.L/N;
		double U[] = new double[N+1];
		
		int i = 0;
		for(double x = data.x0; x <= data.L; x += h)
		{
			if(i == N+1) break;
			
			U[i] = data.uReal(x, T);
			++i;
		}
		
		return U;
	}
	
}
